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A polymer network can imbibe water, forming an aggregate called hydrogel, and undergo 
large and inhomogeneous deformation with external mechanical constraint. Due to the 
large deformation, nonlinearity plays a crucial role, which also causes the mathematical 
difficulty for obtaining analytical solutions. Based on an existing model for equilibrium 
states of a swollen hydrogel with a core-shell structure, this paper seeks analytical so- 
lutions of the deformations by perturbation methods for three cases, i.e. free-swelling, 
nearly free-swelling and general inhomogeneous swelling. Particularly for the general in- 
homogeneous swelling, we introduce an extended method of matched asymptotics to con- 
struct the analytical solution of the governing nonlinear second-order variable-coefficient 
differential equation. The analytical solution captures the boundary layer behavior of the 
deformation. Also, analytical formulas for the radial and hoop stretches and stresses are 
obtained at the two boundary surfaces of the shell, making the influence of the parame- 
ters explicit. An interesting finding is that the deformation is characterized by a single 
material parameter (called the hydrogel deformation constant) , although the free-energy 
function for the hydrogel contains two material parameters. Comparisons with numerical 
solutions are also made and good agreements are found. 

Keywords: Hydrogel; Swelling; Shell; Asymptotic method; Analytical solution; Boundary 
layer. 
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1. Introduction 

Gels, known as a cross-linked solution,!!^ consist of a solid three-dimensional net- 
work of polymer that spans the volume of a liquid medium and imbibes the solvent 
molecules through surface tension effects. When the solvent happens to be water, 
the aggregate is called hydrogel (e.g. edible jelly), which can undergo large and 
reversible volumetric deformation by absorbing or expelling water in response to 
various external stimuli (e.g. temperature, physical or chemical stimuli like light 
and pH). It undergoes a homogeneous deformation without external mechanical 
constraint, but an inhomogeneous and anisotropic one under external constraints 
(often present in practice). 

This paper deals with a core-shell structure, with a shell of gel fixed to a hard 
core of another material (like metal or another polymer), which defines an inner 
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boundary of the network. Due to the good properties such as stabihty, ease of syn- 
thesis, thermalsensitivity and biocompatible nature etc ., such a hydrogel shell has 



such a structure in recent years. It was found that there exists a density fluctua- 
tion within the network which indicates the spatial inhomogeneity. Sometimes, the 
partial detachment of the shell, which means the large stress at the inner surface 
due to the strong swelling, was observed. Thus a good understanding of equilibrium 
swelling states is of crucial importance. However, few analytical results exist for 
such inhomogeneous swelling and consequently there lacks the interpretation of the 
influence of the material parameters on the deformation. 

Equilibrium theories of heterogeneous substances date back to Gibbs,-'^'^ who 
formulated a theory for the inhomogeneous equilibrium state of large deformation 
of an elastic solid in a solvent. Recently extensive studies have concentrated on the 
swelling of gels.ISHSIII] Parti cula rly based on the field theory of Gibbs"'^^ and the 
poroelasticity theory of Biot ,13151 Hong et al.D^ formulated a theory of couple mass 
transport and deformation in gels by considering both the mixing and stretching 
processes, leaving open the free-energy function. For the specific core-shell struc- 
ture of hydrogel, Zhao et al.^^ and Hong et al.^I^ adopted the free-energy function 
introduced by Flory and RehneJl^l and obtained some numerical results of the in- 
homogeneous swelling states, showing large stresses near the core-shell interface. 

The present work is restricted to the equilibrium swelling states (i.e. the long- 
time limit) without considering kinetics, of which the deformation of the network 
is governed by a boundary-value problem. The object of this paper is to seek ana- 
lytical solutions of radially symmetric deformations for such a core-shell structure, 
based on the existing model in Ref. [27l Usually, it is very difficult to obtain ana- 
lytical solutions for an inhomogeneous state of a hydrogel due to the nonlinearity 
caused by the large deformation. In the case of uniform swellings (the water con- 
centratio n is uniform ly distributed), a number of analytical solutions have been 
obtained.^^SElESESl por the present problem, the water concentration is nonuni- 
form, and as far as we know, analytical solutions for this type of problems are not 
available in literature. 

Here we intend to construct the asymptotic solutions for the present problems. 
Identifying a small parameter in the governing equations, we analyze the defor- 
mations by perturbation methods. For the homogeneous deformation, by defining a 
hydrogel deformation constant a we can express the free stretch in a simple formula. 
For the general inhomogeneous deformation, we treat it as a boundary-layer prob- 
lem of a nonlinear second-order variable-coefficient differential equation. It turns 
out that there is a boundary layer near the hard core, but the existing method of 
matched asymptotics does not work for the present problem. Here we introduce an 
extended method of matched asymptotics to construct the analytical solution. More 
specifically, this novel methodology involves the introduction of a transition region 
besides the usual inner and outer regions and using a series solution in this region. 



bioseparation' 
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This paper is arranged as fohows. Section 2 briefly recalls the formulation of Zhao 
et alP3 for the hydrogel in the equilibrium state. We then consider in section 3 the 
free-swelling deformation with no external mechanical constraint, and in section 4 
we discuss a near free-swelling deformation with the fixed hoop stretch at the inner 
surface not far from the free stretch. Section 5 discusses a general inhomogeneous 
deformation without such a restriction on the fixed stretch, where an extended 
method of matched asymptotics is introduced to construct the analytical solution. 
Finally some conclusions are drawn. 



2. Governing Equation 



For the structure of a spherical shell, the spherical symmetric deformation of the 
hydrogel is fully specified by a function r(R). In this section we briefiy recall the 
formulation of Zhao et al.l23 for the hydrogel in the equilibrium state. The field 
equation is 



dSr 



2f^=0, 
R 



(2.1) 



where Sr , sg are the nominal stresses in the radial and circumferential directions 
respectively. 

We adopt the free energy function of the hydrogel first introduced by Flory and 
Rehner (see Ref. I13|14|) and follow the notations in Hong et al."^ 



W(F,C) ==W,{F)+WraiC) 
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where F is the deformation gradient, C is the nominal concentration of water (i.e. 
the number of the water molecules per reference volume in the current state), N 
is the number of polymer chains per reference volume of dry network, kT is the 
temperature in the unit of energy {k is the Boltzmann's constant), Ai,A2 and A3 
are the three principal stretches, and v is the volume per solvent molecule (water), 
X is a parameter from the heat of mixing. The two dimensionless parameters x ^-^d 
vN vary in the ranges 0.1 — 0.5 and 10^^ — 10~^ respectively according to Zhao et 
al.^23 (1/d_/V actually is the number of water molecules occupied the same volume 
of per polymer chain) . 

For a spherically symmetrical deformation, it is easy to deduce from p.2p that 
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(2.3) 
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where \g and are respectively the stretches in the circumferential and radial 
directions, and vC represents the change in volume of the gel, which are given by 

^'^li-' ^'-^ZS' "^ = ^eA.-l. (2.4) 
Substituting (|2. 3112. 41) into (|2.ip . a nonlinear second-order variable-coefficient 

diflterential equation for r{R) arises, which will be solved analytically subjected to 

suitable boundary conditions. 

In the reference configuration (a water-free and stress- free state), suppose that 

the hydrogel shell has the inner and outer radii A and B respectively. Suppose that 

in the current configuration (an equilibrium state immersed in water) the inner 

surface is attached with a rigid core and has the radius r{A) = \qA. At the outer 

surface it is supposed that Sr{B) = or S0{B) = 0. 

Since the change in volume vC is relatively large (see Figure 2(a) in Ref. 27), we 

approximate the term log by the Taylor expansion in terms of ^^^q ■ Then, 

from (P31 - 4Til) we have 
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(2.5) 

+ ■■■ ■ 



We notice that a small parameter vN appears in the equation, so we would like to 
take advantage of this by using perturbation methods to get approximate analytical 
solutions for the following three cases. 



3. Explicit Solution for a Free-svifelling Deformation 

If the hydrogel swells freely with no external mechanical constraint, the deformation 
is homogeneous and isotropic, i.e. Xr = Xe = Xjree = {vCfree + 1)^^'^, which can 
be obtained by solving s,- = (or equivalently sg — 0). Now, we shall deduce the 
explicit asymptotic solution. 

Substituting A := A,. = Ag into (|2.5p we arrive at 



X-X-^^'-^'-^'+...^0. (3.1) 
2vN A4 3vN A7 ^ ' 

Since vCfree is large, we also regard A as a large quantity. From the above equation 

we can see that the term to balance the third term, which is large due to the small 

parameter vN, is the first term A. Thus, they should have the same order, which 

implies that to the leading order 

A=[(l-2x)/(2«7V)]i/5=:«. (3.2) 

We call a to be the hydrogel deformation constant, as we shall see that this single 
parameter plays a dominant role for the deformation. Letting A = aX and seeking 
a perturbation expansion solution of (|3.1|) in the form 

X=l + a-^Xi+a-^X2 + 0{a-^), (3.3) 
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where a is treated as a large parameter, we obtain the formula 

A/^ee = A = a + ia"^ + 0(a"2). (3.4) 
5 

We can see that the single parameter a, which is a combination of the original 
parameters x ^nd vN, determines the deformation (up to the order 0{a~^)), i.e., 
the deformation is not really two-parameter dependent but rather is mainly one- 
parameter dependent. 

Thus the current volume per reference volume is 

1 + vC = = a^ + la. (3.5) 
5 

To the leading order, this result implies that this volume depends on 1/vN by the 
power 3/5, which is consistent with a result obtained before (see eq(13) in Ref. [13]). 
Here, the correction term (|a) is also provided. 

Actually Xfree can be calculated numerically directly from the formula in (j2.3p . 
For several sets of parameters we compare the A free values according to our explicit 
solution and the numerical solution in the following table: 



Table 1. Comparison of explicit solution and the numerical solution for Xfree- 



(vN, x) 


a 


numerical solution 


explicit solution 


error 


(lo-^o.2) 


1.97435 


2.12537 


2.07565 


2.3% 


(10-^,0.2) 


3.12913 


3.21502 


3.19305 


0.68% 


(10-*, 0.2) 


4.95934 


5.00872 


4.99967 


0.18% 


(10-^0.2) 


7.86003 


7.88911 


7.88548 


0.05% 


(| X 10-^,0.3) 


4.95934 


5.01302 


4.99967 


0.27% 



We can see that the very simple formula p.4p for A free agrees with the numerical 
solution very well. As vN or x decreases, a increases, and thus the explicit solution 
becomes more accurate. However, even when is not so small the explicit solution 
gives a very good result already (say, in the case of row one — 0.5065 and the 
error is only 2.3%). This often happens for a perturbation expansion solution: In 
theory one needs that the small parameter tends to zero but in practice the result 
can be valid even when the parameter is not so small. 

The fifth row should be compared with the third row. Although the values for 
vN and x ^re different, the single parameter a has the same value in the two cases. 
It can be seen that the values of A free according to the numerical solution are also 
almost the same. 

4. Analytical Solution for a Near Free-swelling Deformation 

In practice mechanical constraints at the outer and inner surfaces may be present 
and as a result the deformation is inhomogeneous. In this section we consider the 
case that the inner surface R = A has a fixed radial displacement r{A) — A = 
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\q ■ A ~ A (i.e., the stretch Ag — Aq) and the outer surface is still stress-free in 
the radial direction. It is further supposed that |Ao — A^^eel ^ a or |Ao — a| <C a 
for a large a. For this problem, one would expect that the deformation, although 
inhomogeneous, is close to a free-swelling one as Aq is close to A/^ee- Now, we 
proceed to construct the explicit analytical solution. 

For a deformation close to that of a free swelling state, to the leading order, the 
deformation should be given by r{R) = aR. We make the following transformation: 



r(s) — au{s)R, s 



R-A 
B-A' 



(4.1) 



where s is used as the independent variable of u and r in order to convert the 
domain [A, B] to the unit interval [0, 1]. Then, from (|2.5I) and (|2.ip we arrive at 



u'^{s)[u{s) + {s + a)w'(s)]3_ 
4 



(s + a)u"{s) + 



u^{s)[u{s) + {s + a)w'(s)]2 u'^{s)[u{s) + (s + a)u'(s)]3 



u'{sy 



(4.2) 



a 



2u'{s) 



2u'{s) + {s + a)u"{s) 
[u{s) + {s + a)M'(s)]2 u{s)[u{s) + {s + a)u'{s)] 



0{a- 



0, 



where a ^ A/{B — A) is the ratio of the inner radius to the shell thickness. 

At the outer surface s — 1, the outer boundary condition Sr(l) = implies that 



w(l) + (a + l)w'(l)- 



1 



1 



u2(l)[u(l) + (a + l)u'(l)]' 



(4.3) 



u{l) + {a + l)u'{l) 
At the inner surface s — 0, the boundary condition becomes 



(4.4) 



where Aq = Aq — a is regarded as an 0(1) quantity (so that |Ao — a\ <C a). 

Next we seek a regular perturbation expansion solution by considering the pa- 
rameter a to be large. Since to the leading order u{s) should be 1 for a near free- 
swelling deformation, we let 



u{s) = 1 + a ^ui{s) + a ^1*2(5) + 



(4.5) 



At 0(1), equation (|4.2I) and boundary conditions (j4. 31 - 1X41) are automatically sat- 
isfied. At 0(0;^^), we have from equation (j4.2l) that 

[s + a)u'l{s) + Au[{s) = 0. (4.6) 

Solving this equation and further using boundary conditions (|4.3I — fX^ . we obtain 

ui{s) ^ ci{s + a)~^ + C2, (4.7) 
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where 

^ 5(l + a)3 + 4a3' ^ 5(l + a)3+4a3' ^ ' ' 

At 0(a~^), from equation (|4.2p we obtain 

(s + a)M^'(s) + 44(s) = -12c?(s + 0)"^ (4.9) 
Solving this equation and further using boundary conditions (j4.3l - IT^ . we obtain 

U2{s) =di{s + a)-3 + d2~ lcj{s + a)-^, (4.10) 



where 



and 



_ a^l + a)^{5Mi- M2) _ ia^Ah + (1 + afM2 

5(l + a)3 + 4a3 ' 5(l + a)3+4a3 ' ^ ' 



M. = Hc;„-, M...^-^-J^,^«4■ (4.12) 

By transferring back to the original variable R, up to 0(a^^), the solution is 
given by 

- + ill + + + - sHi) ■ 

where r — r/A and R = R/A. We point out that ci, C2, di and (i2 only depend on 
the geometric parameter a and Xq. 

The analytical solution can provide a lot insight information. First, once again 
we can see that the deformation is mainly characterized by the single hydrogel de- 
formation constant a. Next, we shall present the analytical formulas for the physical 
quantities at the inner and outer surfaces. At the inner surface R = A, from the ana- 
lytical solution, the following simple formulas (valid up to 0{1)) can be immediately 
induced: 

2[5(l + a)3-2a3]A5 



" 5(1 -h a)3 4a3 

SB _ 10[(l-Ha)3 -h2a3]A5 Sr _ + af - a^\Xl 



(4.14) 



NkT 5(1 -h a)3 4a3 ' NkT 5(l + a)3-K4a3 

At the outer surface R — B, we have 

9a3AS , 6a3AS 



5(l + a)3+4a3' ' 5(l + a)3 + 4a3' 

sg _ 30a^X* Sr _p ^ ' 



NkT 5(1 + a)3 + 4a3 ' NkT 

The stress values at the inner surface are of particular interest as debonding may 
happen there. We notice that at the inner surface both stress values are proportional 
to the value Aq = Ao — a, the difference between the given stretch and the hydrogel 
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deformation constant a, with the proportional constants dependent on the single 
geometric parameter a, the ratio of the inner radius to the shell thickness. At the 
outer surface, sg is not zero, rather it is an 0(1) quantity proportional to Aq. This 
implies that certain stress in the circumferential direction has to be applied to 
maintain this spherically symmetric deformation. 

To further examine the influence of the geometric parameter a, we consider two 
special situations: a <^ 1 and a ^ 1, which correspond to the cases of the shell 
being very thick and very thin (relative to the inner radius) respectively. 

For a <C 1, at i? = A we have 

Ae = a + AS = Ao, \r^a-2\l, ^ 2X*, « -4AS, (4.16) 



and at R = B we have 



Xe^a + la'Xl A. « a - ^a^A^, ^ « Ga^Xl ^ = 0. (4.17) 



In this case, we see that at the inner surface the magnitude of the stress s,. is twice 
that of Sg and their signs are opposite. Also, sg is very small at the outer surface 
(as a <C 1), which implies that little stress in the circumferential direction needs to 
be applied. 

For a ^ 1, at i? = A we have 

Ae = a + A* = Ao, Xr ^ a ~ -A* w —A* w -—An ~ 0, (4.18) 

° 3 ° NkT 3 NkT 3a ° ' ^ ' 



and ai R = B we have 



Afl w a + An = An, Xr ~ a An, — - — ~ — An, — - — — 0. (4.19) 

° 3 °' NkT 3 °' NkT ^ ' 

In this case, the stresses and stretches at the inner and outer surfaces are approxi- 
mately same, which are somehow expected for a thin shell. In contrast to the first 
case, the stretches Xg and A^ at the outer surface differ from a (or Xfree) by an 0(1) 
quantity, and the stress sg at the outer surface is not small but an 0(1) quantity, 
which means that an 0(1) stress needs to be applied at the outer surface for such 
a deformation. 

The nonlinear second-order variable-coefficient differential equation (|2.ip with 
the boundary conditions r{A) — X^A and Sr{B) — can be solved by using a nu- 
merical method. To examine the validity of our analytical solution obtained above, 
we use a shooting method to get the numerical solution and then compare it with 
the analytical one. In Figure 1, the solution curves according to the two methods 
are plotted. 
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(a) The parameter values are vN = 10~"^ and 
X = 0.2 (a = 4.95934), B = 3A and Aq = 4 
(correspondingly Aq ft: —1). 



X 

7.5r 




(c) The parameter values are vN = lO-"* X 2/3 
and X = 0.3 (a = 4.95934), B = 3A and Aq = 4 
(correspondingly Aq fs! —1). 
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(b) The parameter values are vN = 10— ' and 
X = 0.2 (a = 4.95934), B = 3A and Ao = 3.5 
(correspondingly AJ ~ —1.5). 




Fig. 1. Distributions of the stretches (lines -o- are the analytical solutions and black solid lines are 
the numerical solutions). 

For the chosen geometric parameter in Figure 1, we have a — 0.5. Although the 
value of a is not very large, it can be seen that the analytical solution agrees with 
the numerical one very well. Actually, the maximum relative errorf0 of {Xg, A^) are 
only about (0.2%, 1%), (0.2%, 1.5%), (0.3%, 1.3%) and (0.2%, 1.7%) respectively for 
the four figures. 

We also point out that {vN, x) have different values in Figures (a, b) and Fig- 
ures {c,d), but the a value is the same in all cases. So, the analytical solutions in 
Figures (a, c) and Figures (b, d) are the same respectively. The agreement between 
the analytical solutions and numerical ones show that the deformation is mainly de- 
termined by the single hydrogel deformation constant a, although the free-energy 
function contains two material constants {vN,x)- 

Normally when the outer boundary condition Sr{B) = is used, the stress sg is 
not but an 0(1) quantity (see (|4.15P q). Now, we consider the case that sg{B) = 

''it is defined as the maximum error divided by the maximum value i.e. max \y — y\/ max \y\ 
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instead of Sr{B) = 0. In this case, the boundary condition (I4.3P is replaced by 

"(1) - ^TTTTrTTT^T TwTTvT " ""^^ + = 0. (4.20) 

1*3(1) [m(1) + (a + l)u'(l)] it(l) ^ ^ ^ ^ 

One can proceed to construct the perturbation expansion solution as before. The 
solution expression is still given by (j4.13p but now the expressions for the constants 
are replaced by 

' 5(l + a)3-2a3' ' 5(l + a)3-2a3' ^ ' 

_ a3(l + a)3(5Mi-M2) _ jl + afAh - 2a^ Mi 

"^'^ 5(l + a)3-2a3 ' ^ 5{1 + af - 2a^ ' ^ ' 

where Mi is given by (|4.12p i and 



M,^l+ + + lOc^ (4.23) 



lOCi 5ciC2 

3(TT^ ^ (TT^3 

Simple analytical formulas can also be obtained for the stretches and stresses at the 
inner and outer surfaces, and here we omit the details. 



5. Analytical Solution for a General Inhomogeneous Deformation 

In the previous section, we have assumed that at the inner surface the given stretch 
Ao satisfies the constraint | Aq— a| ^ a. In that case, basically the governing equation 
can be linearized around r — aR so the analytical solution can be obtained by 
solving linear differential equations. Now, we shall proceed to construct the solution 
without the above constraint. Instead, it is supposed that the stretch Aq at ii = A is 
far away from a such that Ao is an 0(1) quantity. For this problem, one cannot avoid 
to deal with some nonlinear second-order variable-coefhcient differential equation(s). 

As mentioned before, for a hydrogel the material constant vN is small, so we 
always take a as a large parameter or a^^ as a small parameter. In general, one 
cannot solve a nonlinear differential equation analytically. However, if a small pa- 
rameter is present in the equation, sometimes one can use singular perturbation 
methods to construct asymptotic solutions. But, for those methods to work, usu- 
ally the equation should become degenerate as the small parameter tends to zero, 
say, it becomes a linear equation or it becomes a first-order equation instead of the 
original second-order equation. For the present problem governed by (|4.2p . we see 
that as tends to zero the leading-order equation is still a complicated nonlin- 
ear second-order variable-coefficient equation. This shows that the existing singular 
perturbation methods do not work for this equation. Here, we introduce a novel 
methodology, which is an extension of the method of matched asymptotics, to con- 
struct the analytical solution. 

We consider the case that the thickness of the shell is relatively large (an ex- 
plicit restriction will be provided later on). We first make some observations on 
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the solution structure. For a thick shell, the boundary condition at the inner sur- 
face should not influence a region some distance away from it (we assume that the 
St. Venant's principle applies). So, there is a region containing the outer bound- 
ary point in which the deformation is near a free-swelling one as the stress-free 
condition can then be satisfied automatically (to the leading order). We call this 
region to be the outer region. When |Ao — a] <C a, we see from Figure 1 that in 
a region near the inner surface the stretches change rapidly. It is reasonable to 
expect that when Aq = 0(1) the stretches also change rapidly in this region. In 
other words, there is a boundary layer region near the inner surface, and we call 
this region to be the inner region. This kii id o f structure can also be seen from the 
numerical solutions obtained in Zhao et aL.^^Ilin the standard technique of matched 
asymptotical^ only an inner region and an outer region exist and the equation in 
the latter region is one-order less than that in the former region. By solving the 
equations in both regions separately and then matching the two solutions together 
to determine the integration constants, the asymptotic solution can be obtained. 
However, for the present problem, the leading-order equation in the outer region, 
which can be obtained by setting = in (|4.2I) . is still a second-order differential 
equation, so one cannot simply match the solutions in the inner region and outer 
region directly. To connect them there should be a third region in between, which 
will be called a transition region. We shall use this transition region to connect the 
outer and inner regions. However, a major difficulty arises: In this region one has to 
deal with the full nonlinear second-order variable-coefficient differential equation, 
which is not solvable analytically! We shall overcome this difficulty by using a series 
expansion for the solution in this region whose interval should be small. The details 
are described below. 

(a) Solution in the outer region 

First we consider the outer solution. The governing equation is still (|4.2I) . and 
the boundary condition (14. 3|) can still be used in the outer region. As mentioned 
before, it is expected that in this region the deformation is near a free-swelling one. 
Therefore, we seek a perturbation expansion solution of the form 



Here, the second-order term is set to be 0{a^^), to be consistent with the governing 
equation (|4.2p and boundary condition (|4.3|) . We substitute this expansion into (|4.2|) 
and ()4.3|) . At 0(1), they are automatically satisfied. At 0{a^^), we find that ui(s) 
satisfies ()4.6p . and the solution expression is 



where C2 and C3 are two integration constants. By further using (j4.3p . we obtain 



u{s) = 1 + a ^Mi(s) + • • • . 



(5.1) 



(5.2) 



ui{s) 



(502-l)(l + a)' 
4(s + a)3 



+ C2. 



(5.3) 
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To sum up, the outer solution is given by 



A a 
where C2 is to be determined. 



(5C2-l)(l + a)' ^ C2(s + q) 
4a(s + a)^ a 



(5.4) 



(b) Solution in the inner region 

Next we consider the inner region (i.e. boundary layer), we should examine the 
full equation (j4.2p . To simplify the equation we introduce the variable r by 



— = f = ar(s) 
A ^ ' 



u{s) 



ar{s) 



s + a 



(5.5) 



Then the equation (j4.2|) becomes 



2(s + a)2 



4(s + a)2 



'(^) 



2a2f'(s) 2a^f{s) 



a^r^(s)r'{s) 



yS + a) {s + a)^ 
{a + s)f"{s) -2f'{s) 2 
{a + s)[f' {s)]'^ f(s) 



4(s + a) 
a^f'^{s)[f'{s)]'^ 

= 0. 



(5.6) 



Suppose that in this region the maximum value f^ax — 0{a^'') {k is to be 
determined) and we write r — a^^r. We note that the value of f at the inner surface 
is 0{a~^) due to the condition f = Aq = 0(1). f,nax should be much larger than this 
value due to the rapid increase of r in the boundary layer region. Thus a restriction 
is fc < 1. To reflect the rapid change of f in the boundary layer, we introduce 
the stretching coordinate X — s/(ae), where ae, the parameter characterizing the 
thickness of the boundary layer, is to be determined. Making this change of variables 
to equation (|5.6p . according to the Van Dyke's principle of least degeneracy,!^ we 
find e — a^^'^^^. Since e should be small, we need fc > 0. And, the equation becomes 



+ 0(e) + O(e') + Oie^a'^-^) = 0, 



(5.7) 



where we denote — *° distinguish from r'(s). Since < fc < 1, 0(a^ ) is 
small. If one uses the Van Dyke's principle of least degeneracy for the O(e^) equation, 
it is required that 0{e'^a^''^^) = O(e^), i.e., k = 6/11. Then, the boundary layer 
thickness parameter ae = aa^^"^^^, which needs to be small, say, ae < 0.15. Thus, 
a restriction is a < 0.15a^^^^^ . 

Multiplying both sides by and integrating once, we obtain (to the leading 
order) 



Ci, 



(5. 



where Ci is the integration constant. This is a first-order differential equation. With 
the boundary condition f(0) = Ao, theoretically there is only one constant Ci to be 
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determined. Actually, the solution fin of the above equation can be represented by 

dy 



X = 



(5.9) 



where Aq := Aoa ^^^^ is a, known constant and / is the root of the cubic algebraic 
equation 



y 



(5.10) 



For the present problem we require f — > with y = f > 0. It is easy to show 
that equation (|5.10p has one and only one positive root, which is given by 



fir, Ci 



^[Vl-C?fVl08-lJ 

Z273 



3iy2[^l-C?yVl08-l]'^^ 
^ 3^[l-^l-CjaVl08]'^' 



373- 



-, Ci < 0, 

Ci=0, (5.11) 
Ci > 0. 



In summary, the inner solution is provided by (j5.9p and (j5.1ip with one constant 
Ci to be determined. 



(c) Solution in the transition region 

Now, we consider the transition region, which is used to connect both the inner 
and outer regions. Since f (or \g) is 0{a^^^^) and 0{a) respectively for the inner 
and outer regions, such a transition region is needed to get the whole solution 
in the whole interval. The independent variable s is in the interval [0,1], and we 
schematically represent the three regions in Figure 2. In this figure, the transition 
region is represented by [sq — A, sq + where sq and A are to be determined. 

transition region 
inner region I outer region 

^ ' ^ 

So - A so + A 1 



Fig. 2. The geometric representation of the three regions. 



In this region, we should use the full equation (|5.6p , and we have (to the leading 
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order) 



2{.3 + a)^ 
4(s + a) 



r"(s) + 
4(.s + a)^ 



2aV'(s) 2a^f(s) 
(s + a) (s + a)2 



(5.12) 



0. 



This is a nonlinear second-order variable-coefficient differential equation, which ap- 
pears to be not solvable analytically! To proceed further, we observe the following: 
The whole interval for s is [0, 1], which is divided into three regions: outer region, 
transition region and inner region. Since usually the outer region is large in a singu- 
lar pert urbation problem (this is also evident from the numerical solutions in Zhao 
et al.l2Il), the transition region should only occupy a small subinterval of [0, 1]. Thus, 
for s in the small subinterval [sq — A, sq -f A] the solution of the above nonlinear 
equation can be expanded as a series (as long as f(s) is sufficiently smooth): 



rtran{s) = tq + ri(s - sq) -f r2(s - so)^ +r3{s- Sq)^ + 



(5.13) 



where ri{i — 0,1,2,3) together with sq need to be determined. Substituting this 
expansion into equation ()5.12p . the left hand side becomes a series of (s — sq)- All 
the coefficients of {s — so)"'{n = 0, 1, 2, 3, • • • ) should be zero. From the coefhcients of 
(s — so)*^ and (s — sq), we can obtain two algebraic relations among the undetermined 
coefficients, which are represented as 



/i(so,ro,ri,r2) = 0, 
f2{so,ro,ri,r2,r3) = 0, 



(5.14) 



where the lengthy expressions of /i and /2 are omitted. To have enough relations 
for the determination of all constants, we need to relate the transition solution to 
the outer and inner solutions. 



(d) Determination of the constants through connection conditions 
We have obtained the solution expressions in the outer region, inner region and 
transition region (see equations (|5.4p . (|5.9|) and (|5.13p ). Each of the outer and inner 
solutions contains one constant and the transition solution contains five constants. 
The subinterval [sq — A, sq -I- A] also needs to be found, so we have another constant 
A to determine. Besides equations (|5.14p i we need another six relations for the 
eight constants Ci, C2, so. A, ri(i — 0,1,2,3), which can be obtained by requiring 
r, r', r" are all continuous at sq — A and sq + A, i.e. 

finis) ~ r tran (s) , r-„ (s) = Tj^jj^ (s) , f "„ (s) = r"^a„ (s) , at s 
rtranis) = r^utis), r't^^,^{s) = r;^t(s), r;;„,(s) = r'^^t{s), at s 

To reduce the above six relations into two relations, by using the solution expression 
(|5.13p we rewrite them as 



= So - A, 
= So + A. 



(5.15) 
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n - 2r2A + 3r3A2 + 0(A3) = f,'„(so - A; Ci)/a, 
2r2 - 6r,A + 0{A') = r^U^o - A; Ci)/a, 

(5.16) 

ro+nA + r2A^+r^A^ + OiA^)=routiso + A;C2)/a, 
n + 2r2A + 3r3A2 + OiA^) = ^^.(so + A; C2)/a, 
2r2 + GrgA + ©(A^) = f;;,(so + A; Ca)/^. 



By some simple manipulations, rj(j = 0, 1, 2, 3) can be eliminated, and as a result 
two equations for the four constants sq. A, Ci, C2 are obtained 



- 3f,„ + Srout + 2A (-3f:,„, + Arf„ + 2Afl,) = 0, 

- 3f + 3f<,„i - 2 A (3f + 2 Ar7„ + Af^, ) = 0. 



where the subscripts "in" and "out" represent the value at sq — A and so — A 
respectively. Another two equations for so, A, Ci,C2 are provided by (|5.14P i 9. By 
the Newton's method, these constants can be easily found. 

To get the solution curve, we take the parameters vN = 10"'* and x = 0.2, 
which yields that a ~ 5. For the geometrical parameter we choose two different 
values B = 3A and AA (i.e., a — 0.5 and 1/3). The stretch at the inner surface 
is chosen to be Ao = 1.077 (which is an 0(1) quantity). For such parameters, by 
solving the system of 4 algebraic equations mentioned above, we find 



a = 0.5 : So = 0.2127, A = 0.06265, Ci = -0.008259, C2 -0.2449, 

(5.18) 

a = 1/3 : So = 0.1670, A = 0.04386, Ci = 0.0008190, C2 = 0.01042. 



For such parameters, we have e = Q!~^°/^* = 0.2332. The parameter ae, a measure 
of the magnitude of the boundary layer thickness has the values about 0.12 and 0.08 
for a = 0.5 and 1/3 respectively, which are consistent with the values of sq — A (0.15 
and 0.12). The subintervals of the transition region [0.150,0.275] and [0.123,0.211] 
(for a — 0.5 and a = 1/3 respectively) are indeed small, as observed before. For 
such small intervals, the series solution (|5.13p should be very accurate. 

Finally we compare our analytical solution with the numerical one obtained by 
a shooting method. The solution curves obtained by two methods for the above 
chosen parameters are plotted in Figure 3. 
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(a) The parameter values are vN = 10""* and 
X = 0.2 (q = 4.95934), B = 3A and Aq = 1.077. 



A 




(c) The parameter values are vN = IQ-* X 2/3 
and X = 0.3 (a = 4.95934), B = 3A and Aq = 
1.077. 




(b) The parameter values are vN = 10-" and 
X = 0.2 {a = 4.95934), B = iA and Aq = 1.077. 



A 




(d) The parameter values are vN = lO""* X 2/3 
and X = 0.3 {a = 4.95934), B = 4A and Aq = 
1.077. 



Fig. 3. Distributions of the stretches (lines -o- are the analytical solutions and black solid lines are 
the numerical solutions). 



It can be seen that the analytical solution agrees well with the numerical so- 
lution. Actually, the maximum relative errors of {Xg,Xr) are about (2.7%, 7.7%), 
(2.9%, 8.0%), (2.8%, 8.1%), (3.0%, 8.4%) respectively for the four figures. Keeping 
in mind that the obtained analytical solution is only valid up to 0(1) and the 0(e) 
(= 0(a-i°/")) terms are omitted, it produces reasonable good results already. 

In Figures 3(a, b) and Figures 3(c, d), the values of {vN, x) are different. However, 
in all cases the a value is same. Since the analytical solution depends only on 
the a value, the analytical curves in Figures 3(a, c) and Figures 3(5, d) are the 
same respectively. The agreement between the analytical and numerical solutions 
once again shows that the deformation is mainly determined by the single material 
constant - the hydrogel deformation constant a. 

Now we shall give a more explicit expression than equation (|5.9|) for the inner 
solution. As it can be seen from (|5.18p that Ci is small, we seek a perturbation 
expansion solution of equation (|5.8p of the form 



= r*iX) + CirliX) + 



(5.19) 
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Substituting the above expansion into equation (|5.8p and using the boundary con- 
dition f(0) = Ao, we obtain 

\ 3/5 



hnis) 



5^ 



18^ 



x + a; 



5^ 



5/3 



V5 / ^r- 

5/3 \ / 5V4 



-2/5' 



-A 



5/3 



(5.20) 



where X = s/{ae) = a^^/^s/a. Then Pin can be immediately recovered by rin{s) 



a' 



5/11, 



((s). In Figure 4, we plot the sohition curves of (|5.9p and (|5.20p . It can be 
seen that the difference is very smalL This supports the validity of the more explicit 
expression (|5.20p . 




(a) The parameter values are vN 
X = 0.2, B = 3A and Ao = 1.077. 




1.05 1.1 1.15 1.2 1.25 1.3 1.35 



R/A 



(b) The parameter values are ,;7V = 10-" and 
X = 0.2, B = 4A and Aq = 1.077. 



Fig. 4. Comparison of the two inner solutions (dots are the solution 115.91 1 and solid lines are the 
explicit solution | |5.20| |). 



Since we have obtained the analytical solution, some simple approximate analyt- 
ical formulas for important physical quantities can be deduced. As in the previous 
section, we consider the stresses and stretches at the inner and outer surfaces. At 
the inner surface R — A, we have 

, _ , , _ ^^5/3 Cia35/33^2/3 
' - - + ^ ' 



+ 



NkT 12 Ao ' NkT 2^A^/^ 2^ ' 

We can see that radial stretch Ar and radial stress Sr/NkT are 0(a^/^) quantities, 
while the circumferential stress sg/NkT is much larger, an 0{a^^^^) quantity. What 
is more, if Ci-terms are neglected (Ci is small), all the stresses and stretches at the 
inner surface are unaffected by the geometric parameter a, the ratio of the inner 
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radius to the thickness. It should be pointed out that, ahhough Ci is small, its 
value depends on a. In the above cases, the Ci -terms in (|5.2ip only have a minor 
influence (less than 0.3%, compared with the flrst terms). 
At the outer surface R = B, we have 

_i 1 - 9C2 , ^ -1 1 - 3C2 
Ao ~ a — a , \r = a + a , 

Sg _-3(l-5C2)_i Sr _p ^^-^^^ 



NkT 2 ' NkT 

Comparing with equation p.4p . we see that the stretch is close to but a little 
smaller than A/^ee while the stretch A^ is close to but a little larger than Xfree- 
The circumferential stress se/{NkT) is of 0(q!~^), which implies that a very small 
stress needs to be applied to maintain this spherically symmetric deformation. 

If we impose the boundary condition se{B) = instead of Sr{B) = 0, the analyt- 
ical solution can be constructed by the same procedure described above. Actually, 
in this case only the expression of the outer solution changes to 

~ .^ s + a _i [ (l-5C2)(l + a)' , ^2(^ + ^) 1 

ToutKs) = a ha — — ■ — -5 \ . (5.23) 

a [ 2a(s -f ay a 

The expressions of the inner and transition solutions and the connection conditions 
are all the same. Also, the stresses and stretches at the inner surface are still given 
by equation (|5.2ip . And, at the outer surface we have 

Xe = a + , A^ = a -I- q;~^(6C2 — 1) 

iVfcT-'^' W = ^^^^^ " • 

We see that in this case an 0{a~^) tensile stress needs to be applied at the outer 
surface to maintain this deformation. 



6. Conclusions 

We study analytically three cases of hydrogel swelling for a core-shell structure, i.e. 
a free-swelling deformation, a near-free swelling deformation and a general inhomo- 
geneous deformation. The hydrogel deformation constant a, which is a combination 
of the two material parameters vN and x, is identified, and it is found that this 
single material parameter plays a dominant role for the deformations in all three 
cases. For the free swelling deformation, a simple formula for the stretch Xfree is 
obtained in terms of a up to 0{a~^). For the near-free swelling one, we obtain the 
analytic solution for the whole region. Some analytical formulas for the stresses and 
stretches at the inner and outer surfaces are given, and it turns out they depend 
linearly on the value Aq — a (where Aq is the given stretch at the inner surface). In 
this case, for a thick shell (the ratio a of the inner radius to the shell thickness is 
small), the geometrical parameter a has little effect. When the shell is thin it is not 
stress-free in the circumferential direction at the outer surface, which indicates the 
boundary conditions s,.(_B) — and S0{B) = are not equivalent. For the general 
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inhomogeneous one, we treat it as a boundary layer problem. An extended method 
o£ matched asymptotics is introduced to solve this problem. More specifically a tran- 
sition region is introduced to connect the inner region and outer region. Further, we 
seek a series solution in the transition region and impose proper connections with 
the inner and outer solutions. Then, theoretically the problem is reduced to solve a 
system of 4 algebra equations. Analytical formulas for the radial and hoop stresses 
and stretches at the inner surface are obtained. It is found that both the radial and 
hoop stresses are very large and the former is an 0(a^/^) quantity while the latter 
is an 0(a^°/^) quantity. Also, these quantities are independent of the geometric 
parameter a (to the leading order). Numerical comparisons are also performed, and 
the results are in good agreement with the analytical ones. 
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